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Unstructured LES of reacting multiphase flows in 
realistic gas turbine combustors 

By Frank Ham, Sourabh Apte, Gianluca Iaccarino, Xiaohua Wu, Marcus 
Herrmann, George Constantinescuf, Krishnan Maheshf AND Parviz Moin 


1. Motivation and objectives 

As part of the Accelerated Strategic Computing Initiative (ASCI) program, an accu- 
rate and robust simulation tool is being developed to perform high-fidelity LES studies 
of multiphase, multiscale turbulent reacting flows in aircraft gas turbine combustor con- 
figurations using hybrid unstructured grids. In the combustor, pressurized gas from the 
upstream compressor is reacted with atomized liquid fuel to produce the combustion 
products that drive the downstream turbine. The Large Eddy Simulation (LES) ap- 
proach is used to simulate the combustor because of its demonstrated superiority over 
RANS in predicting turbulent mixing, which is central to combustion. 

CDP is the flagship LES code being developed by the combustor group to perform LES 
of reacting multiphase flow in complex geometry. CDP is named after the late Charles 
David Pierce (1969-2002) who made several lasting contributions to the LES of react- 
ing flows and to this ASCI program. CDP is a parallel unstructured finite-volume code 
developed specifically for LES of variable density low Mach- number flows. It is written 
in Fortran 90, uses MPI, and has integrated combustion and spray modules. In the first 
five years of ASCI (1997-2002), the principle focus of the combustor group was the de- 
velopment and validation of an entirely new numerical method with the characteristics 
necessary for simultaneously accurate and robust LES on unstructured grids. These com- 
peting ends were both achieved by developing a method around the principle of discrete 
kinetic energy conservation (Mahesh et al. 2003). This numerical algorithm has now been 
extended to variable density, low-Mach number multiphase reacting flows. 

In March of 2003, with the underlying numerics established and validated in a variety 
of flows, a major redesign and rewrite of the code was initiated. Called CDP-a, this 
new version of the code will have all the capabilities of the current code along with 
considerable improvements, all of which are considered critical to achieving the stated 
ASCI goal of a “major advance in engine simulation technology.” 

This paper summarizes the accomplishments of the combustor group over the past 
year, concentrating mainly on the two major milestones achieved this year: 

• Large scale simulation: A major rewrite and redesign of the flagship unstructured 
LES code has allowed the group to perform large eddy simulations of the complete 
combustor geometry (all 18 injectors) with over 100 million control volumes. 

• Multi-physics simulation in complex geometry: The first multi-physics simulations 
including fuel spray breakup, coalescence, evaporation, and combustion are now being 
performed in a single periodic sector (l/18 t/l ) of an actual Pratt & Whitney combustor 
geometry. 
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2. New CDP Code Development: CDP-a 

In March of 2003, development of a new version of CDP, CDP-a, was initiated with 
the goal of performing large-scale simulation of the full Pratt & Whitney combustor. 
When completed, CDP=a will have all the capabilities of the existing CDP code, with 
the following improvements: 

• Parallel preprocessing to allow unstructured simulations of 100 million control vol- 
umes or more, 

• Faster, more scalable solvers based on geometric and/or algebraic multigrid meth- 
ods, 

• Particle/mesh load balancing capabilities, 

• Integrated postprocessing capabilities, including the parallel computation and writ- 
ing of plane and isosurface cuts of instantaneous or statistical data, and the seeding and 
tracking of massless Lagrangian particles, 

• Generally improved code modularity and design, allowing users several levels of 
interaction with the code, and rapid implementation and testing of new models and 
capabilities, 

• Regression testing and version control. 

At the time of writing, the parallel preprocessor has been completed, and a cold flow 
(incompressible) version of CDP-a has been completed and verified against the old code 
in a variety of incompressible flow calculations. A Fortran 90 interface to the algebraic 
multigrid solver of LLNL’s CASC group (Falgout et al. 2002) has been added as an option 
for solving the Poisson system, resulting in a significant overall speedup, although at a 
substantial memory overhead. Particle/mesh load balancing capabilities have been added 
and tested using model particle distributions. Details of these capabilities are described 
in the following subsections. 


2.1. Parallel Preprocessing 

CDP requires as input N p separate grid partition files, where N p is the number of pro- 
cessors that will be used for the LES. These partition files describe the grid (nodes, faces, 
control volumes) associated with a particular partition, and also contain the connectivity 
graph. The preparation of these N p partition files is called preprocessing, and consists of 
four distinct steps: 

• geometry definition 

• grid generation 

• grid partitioning 

• reordering/redistribution 

Up to now, the combustor group simulations have used a single-processor preprocessing 
strategy to generate and preprocess hybrid unstructured grids up to a maximum size 
of about 14 M control volumes (cv’s). This is approximately the maximum size that 
can fit within the 8 GByte physical memory of the high-end workstation used for grid 
generation. At this size, the grid is extremely cumbersome to inspect for quality, and a 
more reasonable limit is about 5 M cv’s. 

To preprocess grids of order 100 M cv’s, it was necessary to develop a new parallel 
preprocessor. Some recent aircraft simulations in the literature have used unstructured 
tetrahedral grids close to this size (25 to 60 M vertices) (Mavriplis & Pirzadeh 1999). 
These grids were generated using a two step approach. First an unstructured grid of 
several million vertices was generated and partitioned on a work-station. These parti- 
tions were subsequently refined using homothetic refinement on a supercomputer. This 
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Figure 1 . Flowchart for the Parallel Preprocessor based on Fluent’s GAMBIT software for 
mesh generation, and our own preprocessing program called PPAREM (Parallel PArtition and 
REorder Mesh). 


approach significantly reduces the cost of the initial grid generation, but gives up in- 
dividual element control. This approach also requires a close coupling to the problem 
geometry during the refinement of elements along the boundary to ensure the grid re- 
mains boundary-conforming. 

The present combustor simulations contain substantially more geometric complexity 
than the aircraft wing-body geometries typically used for simulation. Consequently, en- 
suring boundary conformity during refinement presents a significant challenge. Addition- 
ally, there are many regions of the geometry where control of the final grid distribution 
is desired, such as the geometric details of the injectors or dilution holes. Homothetic 
refinement without subsequent grid smoothing can introduce step-changes in the grid 
spacing that may adversely affect LES accuracy. For these reasons, the present parallel 
preprocessor was developed based on the idea of decomposing the combustor geometry 
into several geometrically-separated regions, and then separately generating the desired 
grid in each of these regions, ensuring only that shared surfaces have identical meshes. 
These meshes are then reassembled on the parallel computer using fast octree searches 
on the common face geometry. 

Figure 1 illustrates the process of mesh generation using the CDP parallel preproces- 
sor. Grid generation is still performed on a single high-end workstation using Fluent’s 
commercial mesh generation software, GAMBIT. For mesh sizes less than about 5 mil- 
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Table 1. Parallel preprocessing of 35 M cv grid 


exact cv count 

35,140,896 

number of msh files, A m 

18 

number of final partitions, N p 

352 

I/O: 


input msh files 

3.72 GB 

output partition files 

3.01 GB 

Wall clock time: 


PPAREM on N m processors: 


read/parse msh files 

504 s 

octree reconnect geometry 

28 s 

N p partition (ParMETIS) 

45 s 

write N p partition files 

1436 s 

Total wall clock time 

2013 s (34 min) 


lion control volumes, the preprocessing is done in an entirely serial way, much as before. 
For larger meshes, however, the geometric decomposition approach described earlier is 
used, and GAMBIT is used to generate a number of separate mesh files. Our parallel 
preprocessing program PPAREM (Parallel PArtition and REorder Mesh) then reads and 
reconnects the global mesh, produces a high-quality partition that is independent of the 
present geometric decomposition using Par METIS, and reorders and cooperatively writes 
the partition files using the MPI-I/O routines. Cooperative writing of the partition files 
represented a challenging programming task, but is the only theoretically scalable way 
to write the final partition. 

The parallel preprocessor was used to prepare two large grids on LLNL’s Frost for the 
cold flow full combustor simulations described in the next subsection. Grid sizes were 
approximately 35 M and 100 M cv’s respectively. Due to the memory requirement of 
ParMETIS (approx. 0.4 GB per million cv’s out of the total requirement of 0.5 GB per 
million CV’s for PPAREM), it was necessary in the case of the 100 M cv grid to produce 
an intermediate partition using a low-memory but poor quality partitioning algorithm 
(simply binning the cv’s based on global index). PPAREM was then run a second time 
on a larger number of processors using this poor quality partition as input. With more 
memory available, a high-quality partition was then calculated using ParMETIS. An al- 
ternative and presumably faster approach would be to increase the memory per processor 
by using more nodes. This was tried unsuccessfully up to 6 nodes for the 100 M case (i.e. 
3 processors per node), at which point the “successive partitioning” strategy described 
previously was adopted. 

Tables 1 and 2 summarize some key indices measured during preprocessing. The most 
expensive and apparently least scalable component of preprocessing is the cooperative 
reordering and writing of the partition files using MPI-I/O. In the worst case, however, 
preprocessing of the 100 M cv grid with 1024 partitions required 4.5 hours. 

2.2. Large-Scale Simulation 

In the process of verification and validation of CDP-a, a variety of incompressible cold 
flow simulations have been performed. The largest of these have been 35 M cv and 100 M 
cv simulations of the full P&W combustor (all 18 injectors). Figure 2 presents a snapshot 
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Table 2. Parallel preprocessing of 100 M cv grid - 2 cases 


exact cv count 
number of msh files, N m 
poor quality partitions, N p i 
final partitions, N p 

I/O: 

input msh files 

poor quality partition files 

final partition files 

Wall clock time: 

1. PPAREM on N m processors: 

read/parse msh files 
octree reconnect geometry 
poor quality partition 
write poor partition files 

2. PPAREM on N. p > processors: 

read poor quality partition 
partition (ParMETIS) 
write N p partition files 
Total wall clock time 


100,241,208 

100,241,208 

18 

18 

96 

96 

512 

1024 

14.28 GB 

14.28 GB 

9.65 GB 

9.65 GB 

8.86 GB 

8.78 GB 

1550 s 

1550 s 

75 s 

75 s 

2 s 

2 s 

1538 s 

1538 s 

122 s 

183 s 

387 s 

501 s 

6445 s 

12385 s 

10121 s (2.8 h) 

16236 s (4.5 h) 



Figure 2. Geometry and contours of the instantaneous velocity magnitude on a plane cut 
through the full P&W 6000 combustor simulation, 35 M cv’s. 


of this simulation illustrating its substantial geometric complexity and fidelity. The max- 
imum ratio of length scales (ratio of overall domain size to the smallest control volume 
sizes) is about 10 4 . When the new code is completed and contains all the combustion and 
spray capability of the existing CDP, the ability to run such large simulations will al- 
low the investigation of important physical phenomena, including azimuthal combustion 
instabilities. Their purpose at this point, however, is mainly for computer science. 
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2.2.1. Multigrid 

Multigrid solvers will play an important role in achieving speed and scalability for the 
very large scale simulations made possible with the new code. As part of the CDP-a 
development, an interface was written to the hypre solvers of LLNL’s CASC group (Fal- 
gout et al. 2002), and specifically the algebraic multigrid solver (AMG). Tests on large 
scale cold-flow simulations (35 M cv’s on between 160 and 512 processors of Frost) where 
the existing preconditioned conjugate gradient solver (PCG) for pressure was replaced 
by AMG yielded a consistent overall speedup of over 4 times. At the largest problem 
size investigated of 100 M cv’s, however, it was not possible to use AMG due to memory 
allocation errors, even when all 64 nodes of Frost were dedicated to this problem . Based 
on these results, the following modifications were made to the solver-related research in 
the combustor group: 

• Although the existing CDP code will be phased out in the next few months as the 
combustion and spray capabilities are implemented and verified in CDP-a, the CDP-to- 
hypre interface was integrated into the old CDP code, where it has resulted in an overall 
speedup of over 2 x for the multiphysics reacting flow simulations. The effect of AMG 
on both the old and new codes is shown in figure 3. 

• The combustor group’s research into geometric multigrid has been given lower prior- 
ity. While the performance of geometric multigrid should theoretically be superior to the 
algebraic method, both methods have the property of linear scalability, and the promising 
results from hypre’s AMG have meant that combustor group computer science resources 
can be concentrated on other non-scalable bottlenecks associated with CDP, including 
parallel preprocessing and particle/mesh load balancing. 

• LLNL’s CASC group has been involved to help solve the memory problems associ- 
ated with the very large scale simulations, and also to provide expertise in optimizing 
the variety of AMG settings for our simulations. 

2.2.2. Scalability 

A scalability study of CDP-a has been performed on Frost. Figure 4 summarizes these 
results for the 35 M cv full combustor cold flow LES, showing acceptable scalability on 
up to 480 processors. 

Although the full combustor simulation was also run with 100 M cv’s on up to 64 nodes 
of Frost, this was only possible with the significantly slower PCG solver for pressure. 
Investigations of the code scalability at these larger problem sizes and over a greater 
range of processors will be made once the memory problems associated with the AMG 
solver have been solved. 

2.2.3. Code Performance 

Table 3 gives some of the key computer science indices for the new CDP-a code. The 
total memory requirement and total file size scale linearly with problem size, so these 
indices are reported per M cv’s. We note that the AMG memory estimate is based on 
a number of smaller simulations made on 32 processors or fewer. Clearly such an aggre- 
gate memory requirement does not give the complete picture for large-scale simulations. 
Assuming linear scaling, the 100 M cv simulation should require 100 M x 3.5 GB = 350 
GB, less than the physical memory available on 32 and certainly 64 nodes of Frostf . As 
discussed earlier, however, attempts to run the 100 M cv simulation with AMG on both 
32 and 64 nodes of Frost resulted in memory allocation errors. 

j Frost has 64 nodes with 16 GB/node and 16 processors/node 
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□ scalars 

| momentum + pressure 

□ sprays 


□ momentum 
| pressure 



wall clock time per time step [seconds] 


wall clock time per time step [seconds] 


(a) old CDP, multiphysics LES 


(b) new CDP-a, cold flow LES 


Figure 3. Speedup of the old and new CDP codes by changing the pressure solver from pre- 
conditioned conjugate gradient (PCG) to hypre’s algebraic multigrid (AMG): a) Old CDP code, 
P&W reacting flow single sector simulation with 2 M cv’s on 80 processors of Frost (« 25, 000 
cv’s/processor). b) New CDP-a, cold flow LES of full combustor, 35 M cv’s on 480 processors 
of Frost (« 73,000 cv’s/processor). 



Processors, N p 

Figure 4. Scalability of the new CDP-a code on Frost (with hypre AMG solver for pressure). 
Problem is the cold flow LES of the full combustor geometry with 35 M cv’s. Scaling is reported 
relative to the 160 processor case, which was the smallest number of Frost processors that could 
run this problem. 


2.3. Particle/Mesh Load Balancing 

The geometric locality of particles combined with the sequential nature of particle/mesh 
time advancement can create load imbalance. By integrating components of the parallel 
preprocessor with CDP-a, it was straight-forward to add particle/mesh load balancing 
capabilities to the new code. Because the code has very long run times relative to the 
startup time (i.e. reading the partition and restart files), we developed the load balancing 
strategy to work like a restart. When a threshold value of load imbalance is reached, a 
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Table 3. Some computer science indices for CDP-a 


Total memory requirement 


- with PCG solver for pressure 

1.5 GB per M cv’s 

- with AMG solver for pressure 

3.5 GB per M cv’s 

I/O (MPI-2): 


- input partition files 

0.10 GB per M cv’s 

- input restart or output result files 

0.20 GB per M cv’s 

Floating point performance (Frost): 


- sustained Megaflops per processor 

78.9 Mflops 

- percentage of peak 

5.3 % 

- flops to memory references 

0.294 

Sample wall clock time (35 M cv on 


480 processors of Frost, cold flow): 


- startup (read partition & restart files) 

13 min 

- 20,000 time steps at 11 s/step 

61 h 

- writing result/restart files 

5 min 

- cooperative writing 2D plane data file 

1 min 


new partition is calculated using the multi-constraint option available in ParMETIS. 
This new partition will have (approximately) the same number of cv’s and particles 
on each processor. This new partition is then cooperatively written to files using the 
routines developed for the parallel preprocessor. The simulation is then stopped, and can 
be restarted using the new partition. 

At the time of this writing, these capabilities have been tested on small problems using 
model particle cost distributions only. Figure 5 illustrates a model particle cost function 
and an initial 32-partitioning produced by the parallel preprocessor using ParMETIS in 
single-constraint mode. ParMETIS was then used to repartition the problem with the 
dual constraints of equal numbers of cv’s and particles on each partition. The theoretical 
simulation timings for the initial partition and multi-constraint partition are shown in fig- 
ure 6 a) and b) respectively. Times are normalized relative to the perfectly balanced case. 
For this particular particle distribution, the initial partition has a normalized computa- 
tion time of over 4. Using a multi-constraint repartitioning, the normalized computation 
time is reduced to 1.04. This implies a speedup of about 4 x, assuming no change in 
parallel efficiency. The multi-constrained partitioning does, however, increase the edge 
cut - in this case by about 40%. The extent to which this reduces the theoretical speedup 
will be determined when this load balancing capability is tested on real problems. 


3. Spray Models in CDP 

The objective of the spray simulation effort is to develop a numerical framework based 
on Lagrangian particle models to perform high-fidelity LES of reacting multiphase flows 
encountered in realistic gas-turbine combustion chambers. Three regimes of flow devel- 
opment are commonly observed inside these combustors (as shown in Fig. 7): 1) ’primary 
breakup regime’ and ‘dense regime’, where the liquid film exhibits large scale coherent 
structures that interact with the gas-phase and disintegrate into filaments by forming 
Kelvin-Helmholtz type instability, 2) ‘intermediate regime’ where the liquid blobs formed 
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(a) model particle cost (b) initial 32-partitioning 

Figure 5. Sample problem used to demonstrate particle/mesh load balancing capabilities in 

CDP-a. 
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time 3 
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(a) single-constraint 

32-partition, edge cut = 
18038, normalized time = 
4.05 


(b) multiple-constraint 32- 

partition, edge cut = 23043, 
normalized time = 1.04 


Figure 6. Computation time on each processor for sample particle/mesh load balancing 

problem shown in figure 5. 


undergo secondary breakup, and 3) ‘dilute regime’ where the droplets evaporate, the fuel 
vapor mixes with the surrounding hot gas giving turbulent spray flames. Droplet defor- 
mation and collision are also important features in the intermediate regime. Figure 7 
summarizes the characterization of the spray formation into the three regimes, based on 
the ratio of the characteristic length scale of the liquid to the available grid resolution 
i/ Ax and the liquid phase volume fraction 0 p . 

The main task was to integrate models capturing these complex phenomena into the 
unstructured LES code (CDP), perform comprehensive validation simulations of each 
model, and initiate large-scale turbulent, multiphase, reacting flow simulation in realistic 
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intermediate regime 
dense regime _ . . 



0 p «l 

Figure 7. Regimes of spray formation defined by the ratio of the characteristic length scale of 
the liquid, £, to the local grid spacing, Aa: and the liquid phase volume fraction, 0 P . 

combustor geometry. A summary of Eulerian/Lagrangian equations with two-way cou- 
pling for interphase mass, momentum, and energy transport is given below. The subgrid 
models for sprays including secondary breakup, droplet deformation and drag, droplet 
evaporation, a hybrid particle-parcel approach for sprays, and accounting for finite-size 
effect of droplets in the dense spray regime are described in brief. Also, two novel ap- 
proaches to more accurately simulate the initial ’primary breakup regime’ and the ’dense 
regime’ are briefly outlined. 


3.1. Gas-Phase Equations 

We solve the variable density, low-Mach number equations with two-way coupling be- 
tween the gas-phase and liquid particles. The formulation is based on the flamelet /progress 
variable approach developed by Pierce & Moin (2001) for LES of non-premixed, turbulent 
combustion. The gas-phase continuity, scalar, and momentum equations are, 


dTgZ 

dt 

dpgC 

dt 

dpgUi 

dt 


+ 


+ 


dpgUj _ d Pg , ~ 

dx dt + 

l/thj l J L 

dJGjZuj d dZ 

dxj dxj P g0lz dxj 
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= d^ {f>9ac d^ ] 
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dxi dxj 
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(3.1) 

(3.2) 

(3.3) 

(3.4) 

(3.5) 


_ 1 ,duj dilj 1 du k 

2 d'U j dui 3 dxk 


(3.6) 


where the unclosed transport terms in the momentum and scalar equations are grouped 
into the residual stress qij, and residual scalar flux qzjandqcj ■ The dynamic Smagorinsky 
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model by Moin et al. (1991) is used to close these subgrid terms as demonstrated by 
Pierce & Moin (2001). For a two-fluid (air + fuel) mixture, one conserved scalar (the 
mixture fraction, Z) and anon-conserved scalar (the progress variable, C ) are solved. The 
gas-phase properties such specific heat, molecular weight, density, viscosity, heat-release, 
and source terms in the progress-variable equation, toe, are obtained from lookup tables 
generated using flamelet theory for non-premixed combustion and are dependent on the 
local values of Z , C, and the subgrid mixture fraction fluctuations, Z . 

With the presence of the liquid phase and inter-phase mass, momentum, and energy 
transport, additional source terms are added into the continuity, scalar transport and 
momentum equations. As the droplets evaporate the heat of vaporization is taken from 
the gas-phase and there is evaporative cooling of the surrounding gas. This gives rise to a 
sink term in the energy equation. By assuming adiabatic walls and unity Lewis number, 
the energy and scalar equations have the same boundary conditions and are linearly 
dependent. Only one scalar equation (for mixture fraction, Z) is solved and other scalars 
including temperature are deduced using flamelet tables. The evaporative cooling effect 
(heat of vaporization) is accounted for in the equation of state during the generation of 
the flamelet tables. The heat content of the liquid fuel is taken into account by computing 
an effective gaseous fuel enthalpy and is used in solving the flamelet equations. 

3.2. Modeling Primary Breakup/Dense Regime 

The initial breakup of the turbulent liquid film into large coherent structures is called 
primary breakup. It is dominated by the interaction of the turbulent liquid film with 
the surrounding turbulent gas-phase giving rise to liquid surface instability waves. These 
interfacial instabilities are important in the overall spray evolution and droplet formation 
process. However, the dynamics of the phase interface is highly complex and as of today 
not well understood. It in fact remains one of the outstanding unresolved problems in the 
area of spray simulation how to model this initial breakup in turbulent environments cor- 
rectly. The majority of the atomization modeling effort is based on Lagrangian particle 
tracking and secondary breakup mechanisms as described later. Implicit in these models 
are the assumptions that a) the characteristic length scale of the liquid, £, is small com- 
pared to the grid size, Ax, and b) the liquid volume fraction © p is small. Although these 
models are successful in predicting secondary breakup in the ‘intermediate’ and ‘dilute’ 
regimes, they are not applicable in the ‘primary breakup’ and ‘dense’ regime, because 
their basic assumptions are invalid. 

Two different but in essence complementary approaches are being pursued within the 
group to correctly predict the ‘primary breakup’ and ‘dense’ regime. On the one hand, 
a novel Eulerian based Large Surface Structure model is being developed that explicitly 
tracks the larger scale dynamics of the phase interface, if Ax > 1. On the other hand, 
Lagrangian spray models are being extended to the ‘dense regime’ by accounting for 
drop/particle sizes that are comparable to but still smaller than the grid size, £/Ax < 1. 

3.2.1. Large Surface Structure Model 

The objective of the Large Surface Structure (LSS) Model is to correctly capture the 
dynamics of the phase interface in the primary breakup regime. To this end, we propose to 
follow in principle the LES approach for turbulent flows. All interface dynamics occurring 
on scales larger than the grid size are explicitly resolved and captured by a Eulerian level 
set approach, whereas interface dynamics occurring on the subgrid scales are described 
by an appropriate subgrid model. Details of this model and some preliminary test cases 
and results are given in a separate paper in this annual research brief (Herrmann 2003) . 
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3.2.2. Modeling Finite-Size Droplets 

In a parallel effort, the Lagrangian framework is extended to account for drop/particle 
sizes that are comparable to but still smaller than the grid size, £/ Ax < 1. The equations 
are based on the original spray formulation developed by Dukowicz (1980) which consists 
of Eulerian fluid and Lagrangian particle calculations, and accounts for the displacement 
of the fluid by the particles as well as the momentum interchange between them. The 
gas-phase equations given above are modified to take into account the volume fraction 
(0 = 1 — 0 p ), where Q p represents the volume occupied by the liquid-phase in a given 
control volume. The gas-phase density, p g and pressure p are then replaced by (p g Q) and 
(p0) in the gas-phase governing equations. 

In brief, the Lagrangian particles are advanced using Lagrangian particle tracking. 
The particle volume fraction Q p is then computed on each grid point and then the 
gas-phase equations are advanced. Use of volume fraction, however, adds to numerical 
complexity and accurate numerical schemes conserving liquid mass are being developed. 
Apte et al. (2003c) summarize the theoretical formulation and present some preliminary 
test cases in particulate flows. 

3.3. Modeling Intermediate Regime 

The particle motion is simulated using the Basset-Boussinesq-Oseen (BBO) equations 
(Crowe et al. 1998). It is assumed that the density of the particle is much larger than 
that of the fluid (p p /p g ~ 10 3 ), particle-size is small compared to the turbulence integral 
length scale, and that the effect of shear on particle motion is negligible. The high value 
of density ratio implies that the Basset force and the added mass term are small and are 
therefore neglected. Under these assumptions, the Lagrangian equations governing the 
droplet motions become 


^ = u — u p)+ (l — A)g (3.8) 

where m p is the mass of the droplet, u p the particle velocity components, u gas-phase 
velocities interpolated to the particle location, p p and p g the particle and gas-phase den- 
sities, and g the gravitational acceleration. The drag force on a solid particle is modeled 
using a drag-coefficient, Cd, 


D 


P solid 


3 r Pg |U 9 ~Upl 

4 d Pp d p 


where Cd is obtained from the nonlinear correlation (Crowe et al. 1998) 


(3.9) 


Q = — (i + aRe b p ) . (3.10) 

Here Re p = d p |u g — u p \/p, g is the particle Reynolds number. The above correlation 
is applicable for Re p < 800. The constants a = 0.15,6 = 0.687 yield the drag within 
5% from the standard drag curve. Modifications to the solid particle drag are applied to 
compute drag on a liquid drop and are given below. 
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3.4. Deformation and Drag Models 
In order to quantify the effect of droplet deformation on drag, Helenbrook & Edwards 
(2002) performed detailed resolved simulations of axisymmetric liquid drops in uniform 
gaseous stream. The simulations were performed using an ftp-finite element method (He- 
lenbrook 2002). An unstructured mesh of triangles which deforms with the interface along 
with a dynamic mesh adaptation algorithm was used allowing higher-order accuracy to 
be obtained even though there is a discontinuity at the interface. This gives results for 
the drag which are accurate to within 1%. Based on their computations for a range of 
density and viscosity ratios, range of Weber (We), Ohnesorge (Oh), and Reynolds num- 
bers (Re), a correlation was developed that provides the amount of deformation in the 
form of ellipticity, E, which is defined as the ratio of the height to width of the drop. 
This is given as, 


E = 1 — O.llWe 0 ' 82 + 0.013. l~^-—Oh~°' 55 We 11 (3.11) 

V Pg Pi 

where pi, p g are the viscosities of the liquid and gas and pi, p g are the densities, re- 
spectively. The non-dimensional Weber and Ohnesorge numbers are defined as, We = 
p g U 2 dp/a and Oh = Hi/ Wpi<jd p , where U is the relative velocity between the gas and 
liquid, d p the diameter of the droplet, and er the surface tension. Accordingly, E < 1 
indicates that the drops have more width than height with deformation in the direction 
perpendicular to the relative velocity. These shapes are called oblate shapes. Similarly, 
E > 1 gives elongation in the direction of the relative velocity giving rise to prolate 
shapes. E = 1 implies spherical shapes. 

The effect of droplet deformation is to change the drag force. This effect is modeled 
by using an effective equatorial droplet diameter, d* = d p E -1 / 3 . The particle Reynolds 
number is also modified, Re* = Re p E -1 / 3 (Helenbrook & Edwards 2002). This is used 
in equations (3.9, 3.10) to obtain the modified drag. In addition the effect of internal 
circulation is modeled by changing the drag on a solid sphere as 


Ddrop 
D solid 


f 2 + 3 pi/p g \ /J 

\3 + 3 pi/p g J 


0.03 (H g /pi)Re° p 65 ) 


(3.12) 


3.5. Stochastic Model for Secondary Breakup 
Performing simulations of primary atomization where one tracks the liquid-gas interface 
in realistic combustor geometries is computationally intensive. Development of numeri- 
cal techniques based on hybrid Eulerian/Lagrangian Level-set /Particle tracking schemes 
to capture the primary atomization process is ongoing and will be implemented into 
CDP. However, the current status is to compute the atomization process using advanced 
stochastic secondary breakup models developed (Apte et al. 2003). Emphasis is placed 
on obtaining the correct spray evolution characteristics such as liquid mass flux, spray 
angle, and droplet size distribution. The liquid film is approximated by large drops with 
size equal to the nozzle diameter and undergoes deformation and breakup. The effect of 
high mass-loading on the gas-phase momentum transport is captured through two-way 
coupling between the two phases. 

In this model, the characteristic radius of droplets is assumed to be a time-dependent 
stochastic variable with a given initial size-distribution. The breakup of parent blobs into 
secondary droplets is viewed as the temporal and spatial evolution of this distribution 
function around the parent-droplet size according to the Fokkcr-Planck (FP) differen- 
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tial equation. This distribution function follows a certain long-time behavior, which is 
characterized by the dominant mechanism of breakup: 

= . ( 3 . 13 ) 

where the breakup frequency (u) and time (<) are introduced. Here, T{x,t) is the dis- 
tribution function for x = log(rj), and rj is the droplet radius. Breakup occurs when 
t > tb rea kup = 1/f- The value of the breakup frequency and the critical radius of breakup 
are obtained by the balance between the aerodynamic and surface tension forces. The sec- 
ondary droplets are sampled from the analytical solution of equation (3.13) corresponding 
to the breakup time-scale. The parameters encountered in the FP equation ((£} and (£ 2 )) 
are computed by relating them to the local Weber number for the parent blob, thereby 
accounting for the capillary forces and turbulent properties. As new droplets are formed, 
parent droplets are destroyed and Lagrangian tracking in the physical space is continued 
till further breakup events. 


3.6. Modeling Dilute Regime 

Typical spray simulations do not resolve the temperature and species gradients around 
each droplet to compute the rate of evaporation. Instead, evaporation rates are estimated 
based on quasi-steady analysis of a single isolated drop in a quiescent environment (Faeth 
1977, Faeth 1983). Multiplicative factors are then applied to consider the convective and 
internal circulation effects. We model the droplet evaporation based on a ‘uniform-state’ 
model. The Lagrangian equations governing particle mass and heat transfer processes 
are well summarized by Oefelein (1997) and are described in brief. 


d 

dt. 


( m p) 


= -rh p 


(3.14) 


17lpCpi ~!ft ^ Tp ^ = - T p ) “ mp/ih v (3.15) 

where A h v is the latent heat of vaporization, m p mass of the droplet, T p temperature of 
the droplet, and C Pl the specific heat of liquid. The diameter of the droplet is obtained 

from its mass, d p = {&m p /TTp p ) 1 ^ . h p is the effective heat-transfer coefficient and is 
defined as 


K = ks ( J) /(T g -T s ) (3.16) 

where k s is the effective conductivity of the surrounding gas at the droplet surface. 
The subscript ‘s’ stands for the surface of the droplet. The solutions to the mass and 
temperature equations for quiescent medium are obtained by defining Spalding mass and 
heat transfer numbers and making use of the Clausius-Clapeyron’s equilibrium vapor- 
pressure relationship (Faeth 1977). In addition, several convective correction factors are 
applied to obtain spray evaporation rates at high Reynolds numbers (Faeth 1977, Faeth 
1983). 


3.7. Hybrid Particle- Parcel Technique for Spray Simulations 
Performing spray breakup computations using Lagrangian tracking of each individual 
droplet gives rise to a large number of droplets (ss 10-30 million) very close to the injec- 
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tor. For LES, simulating all droplet trajectories would be ideal, however, this gives rise 
to severe load-balancing problems on parallel processors. We have developed a hybrid 
particle-parcel scheme to effectively reduce the number of particles tracked. A parcel or 
computational particle represents a group of droplets, N par , with similar characteristics 
(diameter, velocity, temperature). During injection, new particles added to the compu- 
tational domain are pure drops (N par = 1). These drops move downstream and undergo 
breakup according to the secondary breakup model and produce new droplets. 

The basic idea behind the hybrid-approach, is to collect all droplets in a particular 
control volume and group them into bins corresponding to their size and other properties 
such as velocity, temperature etc. The droplets in bins are then used to form a parcel 
by conserving mass, momentum and energy. The properties of the parcel are obtained 
by mass-weighted averaging from individual droplets in the bin. For this procedure, only 
those control volumes are considered for which the number of droplets increases above 
a certain threshold value. The number of parcels created would depend on the number 
of bins and the threshold value used to sample them. The parcel thus created then 
undergoes breakup according to the above stochastic sub-grid model, however, does not 
create new parcels. On the other hand, N par is increased and the diameter is decreased by 
mass-conservation. This strategy reduces the total number of computational particles in 
the domain. In a real simulation with breakup and evaporation models, all the particles 
are clustered in a very small region close to the injector. This may still give rise to load- 
imbalance as only a few processors would solve the Lagrangian particle equations. Domain 
decomposition methods taking into account this load imbalance are being developed and 
will be tested in the present large-scale multi-physics simulation. 

4. Validation Studies using CDP 

A systematic and extensive validation effort of CDP code and the Lagrangian spray 
modules was initiated earlier and several simulations were performed to compare the 
LES predictions of gas and liquid/solid phases with the available experimental data. 
Apte et al. (2003a) first performed an LES of particle-laden swirling, co-annular jet and 
compared the results with the experiments of Sommerfeld & Qiu (1991). This validated 
the Lagrangian particle tracking scheme as well as the accuracy of the numerical method 
in swirling flows. In addition, development & validation of the secondary breakup model 
in simplified combustor geometries was also completed (Apte et al. 2003). This spray 
model validation effort was continued and applied to study model predictions for spray 
evaporation and breakup in complex geometries. 

4.1. Validation of Spray Breakup Model in PW Frontend Geometry 

The stochastic model along with the hybrid particle-parcel approach were used to sim- 
ulate spray evolution from the Pratt and Whitney injector nozzle. Figure 8 shows the 
instantaneous snapshot of the spray field along with the injector geometry. The experi- 
mental data set was obtained by mounting the injector in a cylindrical plenum through 
which gas with prescribed mass-flow rate was injected. The gas goes through the main 
and guide swirler to create a swirling jet into the atmosphere. Liquid film is injected 
through the filmer surface which forms an annular ring. The liquid mass-flow rate corre- 
sponds to certain operating conditions of the gas-turbine engine. For this case, 3.2M grid 
points were used with high resolution near the injector. The domain decomposition is 
based on the optimal performance of the Eulerian gas-phase solver on 96 processors. Due 
to breakup, a large number of droplets are created in the vicinity of the injector. With 




Figure 8. Spray evolution from a realistic gas-turbine injector: scatter plot shows the instan- 
taneous droplet locations in z = 0 plane. Large size droplets are injected from the nozzle wall 
in an annular ring to form a hollow cone spray. 



Figure 9. Comparison of radial variation of liquid axial mass flux at various axial locations, 
P&W RANS with TAB model , present LES with stochastic model, o — o exper- 
imental error bar. 


the hybrid approach, the total number of computational particles tracked at station- 
ary state is around 3.5M, which represents approximately 13M droplets. This includes 
around 150,000 parcels. The load-imbalance due to atomization was found to be sig- 
nificant as only l/3 rd of the processors had more than 10,000 computational particles. 
We are looking into advanced load-balancing techniques to reduce this computational 
overhead due to sprays. CDP-(a) has the capability of dynamic load-balancing based on 
particle imbalance as shown earlier. This will be used to perform these simulations to 
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Figure 10. Instantaneous contours of ISO-propyl mass- fraction superimposed by droplet 

locations in 2 = 0 plane. 

significantly improve the efficiency and speedup. It will also facilitate us to investigate 
advanced collision/coalescence models. 

Figure (9) shows comparisons with the experimental data of radial variation of liq- 
uid mass-flowrates using LES with stochastic model for secondary breakup. Also shown 
are the predictions made by Pratt & Whitney’s RANS-based simulation with a Taylor- 
Analogy Breakup (TAB) model. The LES results are generally in good agreement with 
the experiments. The liquid mass- flowrates basically relate to the spray angle and droplet 
size distributions obtained after breakup. The Sauter Mean Diameter (SMD) are well 
captured by the present LES simulation (within 7% of the experimental data). The di- 
ameter distributions when compared with the experimental data indicates a broad size- 
distribution as opposed to the one predicted by RANS with the TAB model. The size 
distribution also indicates presence of large number of small size droplets which could 
be attributed to the lack of models addressing droplet collision and coalescence. In ad- 
dition, the initial droplet size at the injector nozzle is assumed to be constant in these 
simulations, whereas it may vary depending on the local conditions governing primary 
atomization. A further investigation with inclusion of collision models as well as using a 
size distribution at the inlet should be performed in order to accurately predict the spray 
characteristics. 


In order to validate the evaporation model and the variable density formulation in CDP, 
simulation of a coaxial non-swirling jet was performed following the configuration used 
in the experiments by Sommerfeld & Qiu (1998). This flow configuration was chosen 
because of its direct relevance in gas-turbine combustion chambers. In addition, in these 
experiments the boundary conditions for the liquid phase specifying the inlet droplet size 
distribution and their correlation with droplet velocity is well-defined. The gas-phase 
temperatures are not high enough to produce spray flames. This isolates the droplet 
evaporation problem from spray breakup and combustion and is very useful in validating 
the evaporation models used in CDP. Accordingly, only one scalar equation (mixture 
fraction) is solved and the gas-phase temperature is obtained from the ideal gas law. 

The selection of droplet size and velocity at the injection surface is based on the given 
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Figure 11. Comparison of radial variation of evaporation statistics at various axial locations 

with the experimental data of Sommerfeld & Qiu (1998), LES, o experimental data: a) 

mean and rms of axial particle velocity, b) mean axial liquid mass flux, c) mean and rms of 
particle droplet diameter. 


experimental PDF of size- velocity correlations. The grid used consists of 1.5M hexahedral 
cells and around 0.5M particles were present in the computational domain at stationary 
state. Figure (10) shows a snapshot of mixture fraction contours superimposed by scatter 
plot of ISO-propyl alcohol particles in a coaxial combustor. The droplets are injected 
near the inlet circular wall of cross-sectional diameter 40 mm. The droplet velocity- 
size correlation depicts a conical spray with a spray angle of around 60°. Figure (11) 
shows the droplet statistics. The radial variations of mean and rms droplet velocity, 
mean axial mass flux, and mean and rms droplet diameter compared with experimental 
data of Sommerfeld & Qiu (1998) are shown. The profiles are in good agreement with 
the experimental data. The droplet mean diameter profile is typical of the hollow cone 
atomizer, where smaller droplets are found in the core region and larger droplets near the 
edge of the spray. Away from the inlet section, the droplet mean diameter distribution 
becomes more uniform over the cross-section and slowly decreases in the downstream 
direction because of evaporation. The axial mass-flux also decreases towards the exit due 
to evaporation and is well captured by the present simulation. At x = 0 and x = 0.786 
the profiles of droplet mass flux show two-peaks associated to hollow-cone spray. Due to 
the recirculation region downstream of the center-body, the droplet mass-flux becomes 
negative. Further downstream spreading of the spray is hindered by the annular air-jets 
and maximum of the mass flux moves towards the centerline. These features are well 
captured by the present LES computation. 
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5. Large Scale Multi-physics Simulation using CDP 

A multi-scale, multi-physics simulation of turbulent reacting flow in a realistic Pratt 
& Whitney combustor was initiated. This simulation includes all the complex models 
for spray breakup, evaporation, and turbulent combustion described in section(3). This 
simulation will serve as the first validation study of the reacting multiphase flow solver 
(CDP) in complex combustor geometry. Figure (12) shows the combustor geometry in the 
symmetry plane z = 0. The simulation is performed for a single injector which represents 
one sector of the full combustor containing 18 injectors. Liquid fuel (Jet-A) enters the 
combustion chamber through an annular ring at the injector exit. This liquid film is 
approximated by large drops of the size of the annulus radius. These drops are converted 
by the surrounding hot air, they break, evaporate, and the fuel vapor thus formed mixes 
with the surrounding air giving a non-premixed spray flame. The computational grid 
consists of 1.9M hybrid elements (hexes, pyramids, and tets) with fine resolution close 
to the injector. 

The flamelet library for Jet-A fuel at gas-turbine engine operating conditions, was 
generated by using a surrogate fuel of 80% n-Decane and 20% 1-2-4 tri-methyl-benzene 
which closely follows the chemical kinetics and reaction rates of the Jet-A fuel. Around 
1000 elementary reactions among 100 chemical species were used to generate these tables. 
The chemical kinetics of surrogate fuel was compared with original fuel in terms of 
prediction of pollutants in laminar flames and showed good agreement. In this multiphase 
simulation, the flamelet tables are generated by considering the heat of vaporization of 
the liquid fuel. This is obtained by using heat-content for an equivalent gaseous fuel 
which gives large densities in the pure fuel- vapor region ( Z = 1). 

Figure (12) shows instantaneous snapshots of progress variable (C), mixture fraction 
( Z ), normalized temperature (T/To), and velocity magnitude (Umag) in the z = 0 sym- 
metry plane. Scatter plot of droplets forming a conical spray is also shown. A highly 
complex unsteady flame is observed in this simulation. High temperatures in the com- 
bustion chamber are reduced by large mass-influx through the dilution holes. This reduces 
the exit temperature considerably. Strong recirculation zone with conical spray flame is 
observed near the nozzle. The droplets evaporate completely close to injector and do not 
appear beyond ( x = 2). The experimental data available for validation includes pres- 
sure drops across different components, mass-splits through the inner and outer dilution 
holes, and swirlers. The comparison of these quantities with the experimental data was 
shown to be within 10% for the cold flow simulation. Similar results are obtained for 
the reacting flow case. In addition, the exit plane temperature was measured at 5 loca- 
tions. The temperature field is within 10-15% of the experimental data near the center. 
The predicted temperature near the walls is much lower and is still evolving with time 
indicating that the flow has not reached the statistically stationary state. 

This computation was performed on 80 processors on Frost and involved around 0.5M 
computational particles. The algebraic multigrid solver (AMG) developed at Lawrence 
Livermore is used to solve the Poisson equation. The convergence rate of the AMG solver 
was compared with the conjugate gradient solver (PCG). For each time-step, the conju- 
gate gradient solver typically requires 2000 iterations, whereas the AMG solver requires 
10-12 cycles corresponding to 500-600 conjugate gradient iterations. The AMG-solver 
gave an overall speed-up of 2-3 times over the PCG solver and reduced the computa- 
tional time substantially. 
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6. Summary 

CDP is the flagship LES code developed by the combustor group to perform LES 
of reacting multiphase flow in complex geometry. CDP is named after the late Charles 
David Pierce (1969-2002) who made several lasting contributions to the LES of reacting 
flows and to the DOE’s ASCI program. CDP is a parallel unstructured finite- volume code 
developed specifically for LES of variable density low Mach-number flows. It is written in 
Fortran 90, uses MPI, and has integrated combustion and spray modules. A wide range 
of simulations of multiphase, reacting turbulent flows have been performed using CDP. It 
has been shown that the predictive capability of CDP for complex flows in realistic gas- 
turbine configurations is good. A major redesign and rewrite of CDP has been initiated 
to improve and add to the current capabilities. In the past year, major accomplishments 
included the following milestones: 

• Multi-physics simulation in complex geometry: The first multi-physics simulation in- 
cluding fuel spray breakup, coalescence, evaporation, and combustion is being performed 
in a single periodic sector (l/18 t,t ) of an actual Pratt & Whitney combustor geometry. 

• Large scale simulation: Performed large eddy simulations of the complete combustor 
geometry (all 18 injectors) with over 100 million control volumes using CDP-a. 

Further development of CDP-a will include the implementation of advanced com- 
bustion and primary atomization models, and modification of the numerical algorithm 
to capture acoustic waves and combustion instabilities based on the compressible for- 
mulation of Wall et al. (2002). Several large-scale, multiphysics simulations in the PW 
combustor using one or more sectors will be performed to assess the accuracy of the 
overall formulation. 
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Figure 12. Instantaneous snapshots (from top to bottom) of progress variable, C, mixture 
fraction, Z, temperature normalized by inlet air temperature, T/To, and velocity magnitude, 
\/m 2 + v 2 + w 2 superimposed with droplet scatter plot in Pratt & Whitney combustor, z = 0 
plane. The geometry has been distorted. 


